function Im=newBeamForming(eco,pulso,pitch,c,retardo,fs,ResX,ResZ,x0,x1,z0,z1)

%INPUTS
%eco: matriz M2xN, siendo N los canales y M2 el numero de muestras
%pulso: idem M1xN, tiene que estar muestreado a la misma tasa que el eco
%xx, zz: vectores fila que marcan la grilla sobre la que se trabajara
%pitch: distancia entre sensores
%c: vel de propagacion
%fs: frec. de muestreo
%retardo: demora entre comienzo de pulso en emisión y comienzo de recepción

[M2,N]=size(eco);
[M1,N1]=size(pulso);

xx=x0+[0:(ResX-1)]*(x1-x0)/(ResX-1);
zz=z0+[0:(ResZ-1)]*(z1-z0)/(ResZ-1);


%Chequear que las tasas de muestreo sean iguales
%Chequear que el #sensores sean iguales (no necesario)
%Chequear que xx y zz sean del mismo tamaño
%Si son matrices, pasarlos a vectores
xx2=reshape(xx,1,size(xx,1)*size(xx,2));
zz2=reshape(zz,1,size(xx,1)*size(xx,2));
R=size(xx,2);
%chequeo de congruencia entre las profundidades exploradas y las
%disponibles

%extender eco para considerar retardo (agregar algo para evitar aliasing?)
eco_ext=[zeros(ceil(retardo*fs)+M1,N);hilbert(eco)]; %OJO con el M1, depende de si el retardo se cuenta desde el comienzo o el final del pulso. No hay mejor manera de hacerlo?
%extender pulsos al MISMO tamaño y analitizar
newM=size(eco_ext,1);
pulso_ext=zeros(newM,N);
pulso_ext(1:M1,:)=hilbert(pulso);

%Calculo matriz G (newMxN):
G=fft(pulso_ext,[],1);
G=G(1:(newM/2+1),:);
clear pulso_ext;

%Calculo matriz F (newMxR):
F=zeros((newM/2+1),R);
for i=1:N
    aux=calcularHi(i,pitch,xx,zz,c,newM,fs).*(G(:,i)*ones(1,R));
    F=F+aux;
end
clear aux;
clear G;

%Calculo matriz M:
%Mat=zeros((newM/2+1),R,N);
newMat=zeros((newM/2+1)*N,R);
for i=1:N
    %Mat(:,:,i)=calcularHi(i,pitch,xx,zz,c,newM,fs).*F;
    newMat((i-1)*(newM/2+1)+(1:(newM/2+1)),:)=calcularHi(i,pitch,xx,zz,c,newM,fs).*F;
end
clear F;

%Ahora, multiplicando por s en la 2da dimensión pasaría a tener una matriz
%eco esperado newMxN, y quiero minimizar la diferencia entre dicha matriz y
%el eco que tengo.

%Opción 1: concateno en la primera dimensión para todos los i, de manera de
%tener un sistema E=Ms, donde e es vector y M matriz de 2 dim
%newMat=zeros((newM/2+1)*N,R);
E=zeros((newM/2+1)*N,1);
for i=1:N
   %newMat((i-1)*(newM/2+1)+(1:(newM/2+1)),:)=Mat(:,:,i);
   aux=fft(eco_ext,[],1);
   E((i-1)*(newM/2+1)+(1:(newM/2+1)))=aux(1:(newM/2+1));
end
clear Mat;
clear aux;
clear eco_ext;

s=newMat\E; %Resolución a lo groncho!
Im=s;

end

function Hi=calcularHi(canal,pitch,xx,zz,c,M,fs)
%INPUTS
%canal: indica transductor para el que se calcula la transferencia
%pitch: distancia entre transductores
%xx,zz: grilla para la que se calcula, deben ser VECTORES fila del mismo tamaño
%c: vel. de propagacion
%M: numero de muestras totales en frecuencia
%de frecuencias.
%fs: tasa de muestreo

dist=sqrt((xx-(canal-1)*pitch).^2+zz.^2);
delays=fs*dist/c; %En muestras
aux=exp(-2*pi*1i*mod(delays/M,1)); %fila
Hi=matrizPotencia(aux,[0:(M/2)]');
gain=ones(M/2+1,1)*(1./dist);
Hi=Hi.*gain;
end